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Abstract 

Planets are formed from collisional growth of small bodies in a protoplanetary disk. Bodies much 
larger than approximately 1 m are mainly controlled by the gravity of the host star and experience 
weak gas drag; their orbits are mainly expressed by orbital elements: semimajor axes a, 
eccentricities e, and inclinations i, which are modulated by gas drag. In a previous study, d, e, and i 
were analytically derived for e <C 1 and i H/a, where H is the scale height of the disk. Their 
formulae are valid in the early stage of planet formation. However, once massive planets are formed, 
e and i increase greatly. Indeed, some small bodies in the solar system have very large e and i. 
Therefore, in this paper, I analytically derive formulae for d, e, and i for 1 — <C 1 and i H/a 
and fox i ^ H/a. The formulae combined from these limited equations will represent the results of 
orbital integration unless e> 1 or i>7r — if/a. Since the derived formulae are applicable for 
bodies not only in a protoplanetary disk but also in a circumplanetary disk, I discuss the possibility of 
the capture of satellites in a circumplanetary disk using the formulae. 
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Correspondence/Findings 

Introduction 


Planets are formed in a circumstellar disk composed of gas and solid materials (solids are of the order of 
1% in mass). The solid material is initially sub-micron grains, which are controlled by an aero dynamical 
frictional force that is much stronger than the gravity of the central star ( Adachi et al. . 19761 hereafter 
AHN). As solid bodies grow, gas drag becomes relatively less important. Once bodies get much larger 
than 1 m, they have Keplerian orbits around the central star that are slightly altered by gas drag; then, 
their orbits are characterized by orbital elements: semimajor axes a, eccentricities e, and inclinations i. 
These bodi es grow via collisio ns, and the collisional rates are given by relative velocities determined by e 


and i (e.g.. llnaba et al. LuPOlh . Damping due to gas drag and stirring by the largest body in each annulus 


of the disk mainly control e and i, which evolve in the protoplanetary disk during planet formation. In 
addition, radial drift d ue to gas drag, which is expre ssed by d, reduces small bodies, which stalls the 
growth of bodies (e.g., Kobavashi et al. . 2011 . 2010l) . Therefore, the time derivative of a, e, and i (a, e, 
and i) caused by gas drag are very important for planet formation. 


Protoplanets are formed out of collisions with kilometer-sized or larger bodies called planetesimals. 
While protoplanets grow, e and i of planetesimals are determined by the Hill radius of the protoplanets, 
and their e an d i are smaller than 0.3 unless the protoplanets are greater than ten Earth masses (see 
equation 15 of iKobavashi et al. Il2010h . Therefore, AHN derived formulae of d, e, and i due to gas drag 
























for a body with low e ^ 0.3 and i <C 0.1. However, e and i may possibly increase when more massive 
planets are formed. In deed, in the solar system , some comets, asteroids, and Kuiper belt objects have 


very high e and i (e.g.. lKobavashi et al. Ll2005n . In addition, if inclined and eccentric orbits of irregular 


satellites a round Jovian plane ts are originated from captures due to interaction with circumplanetary 


disks (e.g.. lFuiita et al. Ll2013h . these captured bodies with high e and i evolve their orbits in the disks. 


Therefore, analytic formulae for d, e, and i for bodies with high e and i are helpful for the analysis of 
small bodies in the late stage of planet formation. 

In this paper, I first introduce a model f or gaseous disks such as protoplanetary and circumplanetary 
disks, and then, I revisit the derivation of Adachi et al. ( 1976h for the analytic formulae of d, e, and i. 
Next, I derive d, e, and i for bodies with high e and/or high i. By combining these limited solutions, 
I construct approximate formulae for d, e, and i, which are applicable for all e and i unless e > 1 or 
i > TT — H/a. Lastly, I discuss the orbital evolution of satellites captured by circumplanetary disks using 
the derived analytic formulae for d, e, and i. 


Nebula disk model and gas drag law 

In order to evaluate the drag force due to nebula gas, the disk model is set as follows. A gaseous 
disk rotates around a central object with mass M*, which is axisymmetric and in a steady state. In a 
cylindrical coordinate system (r, 9, z), the gas density p is defined from fhe force equilibrium in fhe z 
direction in a vertical isofhermal disk as 
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P = 




■ exp 
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( 1 ) 


where a{r){= pdz) is fhe surface densify of fhe nebula disk, H{r) = is the disk scale 

height, Ok = is the Keplerian angular velocity, and G is the gravitational constant. For 

simplicity, the r-dependences of a and c are assumed as a oc r““, c oc r~^, respectivel y. This relations 
give p oc , where a' = a — /3 + 3/2. In the minimum-mass solar nebula model (IHavashi et al. 


1985h . for example, a = 3/2 an d 0 = 1/4. The ang ular gas velocity Hgas is obtained from the force 


equilibrium in the r direction as (ITanaka et al. Ll2002h 
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In Equation Q, the terms of 0{z^/r^) and higher are ignored. This treatment is valid even for inves¬ 
tigation of the gas drag effect on highly inclined orbits because the gas drag (and the nebula gas) is 
negligible at a high altitude {z^ H). 


At the midplane of the disk, the relative velocity difference between the gas motion and the Keplerian 
rotation is given by 


pir) = 


H JC fin 


Hr 


z=0 


If ^3 

--U + P + - 


H{ry 


( 3 ) 


For a body with mass m and radius d, gas drag force per unit mass can be written as AHN 

Fd = CBTTd‘^puu/2m = Apuu, (4) 

where Cd is the dimensionless gas drag coefficient, u is the relative velocity vector between the body 
and the gas, u =| u |, and A = Gd-kP?/ 2m. Although Cd depends on Mach number M and Reynolds 
number Re, Cd is a constant for high Re{d 1 km in the minimum-mass solar nebula) or for M 3> 1 
(e or i is much larger than H/a) (AHN). 




























General expressions for the change in a, e, and i 


In this paper, I investigate the time variations of semimajor axis o, eccentricity e, and inclination i of a 
body due to gas drag for constant Cd (and then constant A). The time derivatives of a, e, and i are given 
by AHN as 
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where / and uj are the true anomaly and the argument of pericenter, respectively, k = 
re = [1 — sin^(/ + uj) sin^ 
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Pq is the midplane density at r = a, and = (GM^/a)^/^ is the Keplerian velocity. If the variation 
timescales of a, e, and i are much longer than the orbital time, the evolution of a, e, and i follows the 
averaged rate. The orbital averaging is taken as 


\dt) Tk Jo dt ^ 27r/o dt (1 + ecos/)2 


( 10 ) 


where Tk = 27r(a^/GM*)^/^ is the Keplerian period. The same averaging is taken for e and i. 


Case of low e, i, and rj 


For e, t <C 1, lAdachi et al. (1197611 derived the averaged changes in a, e, and i for three cases, (i) p 
e,i, (ii) i ^ e,p, and (iii) i,p, and summed up the leading terms for these cases. This method was 
used to treat u in Equations ([5ll to (|7]l analytically: The assumptions simplify as u ^ p + (e/2) cos / in 
case (i), u « i | cos(/ + cu) | in case (ii), and u k, ey^l — (3/4) cos^ / + {p/2) cos /y^l — (3/4) cos^ / 
in case (iii). Other terms are also simplified, such as p = po(l + a'ecos /). Then, the terms are easily 
averaged over the orbital period by Equation (fTOb . 

The derived formulae are in good agreem ent with the results of orbital integrations for e <C 1 and 
i <C H{a)‘^/a?. While Adachi et al. ( 1976h provided the term of in d, they did not take into account 
the vertical dependence of p, which include s other terms. Si nce the sum of these terms is negligible, 
I thus exclude the term derived by AHN. Ilnaba et al. I(l200lh found that the mean root squares of these 
limited solutions are in better agreement with the results of orbital evolution than the simple summation 




























by Adachi et al. ( 1976h . The averaged variation rates of a, e, and i are therefore given by 
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where K = 2.157 and £' = 1.211 are the first and second complete elliptic integrals of argument y^3/4, 
respectively, and tq = {ApoVY^)~^ is the st opping time due to ga s drag for u = uk- Note tha t I corrected 
an erro r in the factor of the rf' term for e in Adachi et al. (197a), which was pointed out by Karv et afl 
(Il993h . 


For i = 0.01, Equations (fTTl) to (fT3]) are compared with the results of orbital integrations in Figure [T] 
These formulae are valid unless e > 0.2. Moreover, the i dependence in these formulae are valid for 
i < H{a)/a (see Figure |2l). 


Case of high eccentricity and low inclination 


Here, let us consider the case where e is almost equal to unity and i is much smaller than H{a)/a. 
Expanding Equations (l5]l to ([7]) with respect to (1 — e^) under the assumption of / <C H{a)/a, keeping 
only the lowest-order terms of (1 — e^), and applying the orbital averaging such as Equation (fTOl) to 
these equations. 
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The dependences of d and e on / are seen in the integral 'F, while a term proportional to sin 2/ sin 2a; 
in i vanishes by the orbital averaging because of an odd function of /. The integrals of T;, $i, and <1>2 
are functions of a — p. In the minimum-mass solar nebular model, a — /3 is 5/4, and then, T' = 0.79, 
$1 = 0.71, and 4>2 = —0.16. 


The 6 dependences in these formulae are applicable for e > 0.9 as shown in Figure [T] Although the 
effective range of these formulae is limited, the e dependences improve the high e parts in Equations (fTTl) 
to (fT3]) as shown below. 
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Figure 1 The variation rates of a, e, and r as a function of e for r = 0.01 and a; = 7r/2 in the disk with 
H{a)la = 0.1, a = 1.5, and /3 = 0.25. Analytic formulae for low e (gray dotted curves), given by 
Equations (fTTl) to (fT^ . and ones for high e (gray dashed curves), given by Equations (fTdl) to (fl^ . are 
in good agreement with the results of orbital integration (open circles) for low e or high e, respectively. 
The combined formulae (solid curves), given by Equations (1^ to (l3^. are valid for the whole region. 


Case of high inclination 

Next, let us consider highly inclined orbits where ai/H{a) is much larger than unity. Bodies with such 
a high inclination penetrate the nebula disk twice around the ascending and descending nodes through 
an orbital period. Gas drag is effective only around the nodes. Since the body experiences significant 
gas drag around the ascending node (| / + w |<C 1), the leading terms of | / + w | for Equations (l5]l to 
O are 
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Figure 2 The change rates in a, e, and i as a function of i for e = 0.1, and a; = 7r/2 in the same disk 
as Figure [T] Analytic formulae for low i (gray dotted curves), given by Equations (fTTl) to (fT^ . and 
those for high i (gray dashed curves), given by Equations (1^ to (1^ . are in good agreement with the 
results of orbital integration (open circles) for i <C H/a and i S> H/a, respectively. The combined 
formulae (solid curves), given by Equations (1^ to (IMl). represent within a factor of 1.5. 
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Eor this derivation, flgas = since the relative velocity is mainly determined by inclination. 

















In order to apply averaging over half an orbit around the ascending node, a, e, and i are integrated from 
/ = —uj — 7r/2 to / = —cj + 7r/2. Since a, e, and i are Gaussian functions as shown in Equations (l20l) to 
(l22]) . they are negligible for large (/ + and the integral is thus approximated to be that over interval 
[— 00 , 00 ] as follows: 
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where H = H{a). Using this, Equations (l20b to (1221) are integrated around the ascending node, which 
results in the averaged variation rates of a, e, and i in half an orbit. 

The variation rates due to the penetration near the descending node (/ Ri —oj — vr) are obtained in the 
same way as above. Summing up the changes at two penetrations, the averaged changes are given by 
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The validity of Equations (l30l) to (l32l) is shown in Eigures [2] and [3] These formulae are applicable for 
i > 2H/a. 


Combined equations 


The variation rates of a, e, and i in two limited cases for i <C H/a are derived above. The formulae 
for low e do not well reproduce the variation rate in e ~ 1, while high-e formulae overestimate the 
values for low e. Combination of low-eccentricity formulae of Equations (fTTl) to (fT^ with the 1 — 


dependence derived in Equations (fldl) to (fT6l) gives 
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These formulae are given in a very simple way, but they are in good agreement with the results of orbital 
integration if i < H/2a (see Eigures [T]to|3l). 


If H/2a < i < H/a, the variation rates of a, e, and i are obtained from combination of the low-i 
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Figure 3 Same as Figured but for e = 0.9 and dotted lines given by Equations (fT4b to (fTbl) . 


formulae of Equations (1^ to (1^ and the high-i formulae of Equations (1301) to (l3^ . 



where MIN(E), E) is the smaller of D and E. 

In conclusion, the variation rates for a, e, and i are approximately given by 

• Equations (1^ to (1^ for i < El/2a, 

• Equations (1^ to (1^ for the intermediate inclination {H/2a < i < 2H/a), 

• Equations (l30l) to (l3^ for i > 2H/a. 


(36) 
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In the intermediate i, the formulae tend to deviate from the right values but the accuracies are within a 
factor of 1.5 (see Figures [T] to O. It should be noted that these formulae are not applicable to the case of 
i > TT — H/a where a body experiences gas drag with relative velocity ~ vk not only around the nodes 
but also for a whole orbit. 


Application to captured satellites 

Jovian planets have many satellites, which may be formed in circumplanetary disks. Satellites close 
to planets mainly have circular and coplanar orbits and may be formed in the disks. However, distant 
satellites tend to have inclined orbits. Here, I discuss the possibility of the capture of satellites in the 
disks because the formulae for a, e, and i that I derive in this paper are applicable to bodies with high e 
and i. 


Orbital evolution of bodies with high e is predicted from these analytic formulae. When a body is 
captured by gas drag in a circumplanetary disk, e of the captured body is approximately 1. For e > 0.9, 
|e|/e and |d|/a are very large. Variation rate of the pericenter distance q is much smaller than those of a 
and e. Indeed, q = {1 — e)a — ae is estimated to be zero in Equations (fTdl) and (fTSl) . The result is caused 
by the neglect of the higher terms of (1 — e^), and these higher (1 — e^) terms give q/q a positive value 
but q/q is much smaller than |d|/a and |e|/e. Therefore, the orbital evolution occurs along with almost 
constant q. With decreasing e, the orbital evolution changes. Since |d|/a becomes smaller than |e|/e for 
e < 0.5 to 0.6, e decreases with almost constant a. Once e r], a becomes dominant for the orbital 
evolution; the body drifts to the host planet in the timescale of tq/It]^ . 


The bodies that will be sate llites are temporally captured by a planet at first (iSuetsugu and Ohtsuki 


2013fc ISuetsugu et al. 1201 ih . and the apocenter distances of the bod ies decrease to less than the Hill 


radius of the host planet during the temporal capture of bodies (e.g., Fuiita et aL~l 2013 1. The change 
of orbital eccentricity in an orbit around the host planet is given by Ae (e)TK. The body is fully 
captured by gas drag if /capAe ~ 1 during the temporal capture, where /cap is the number of close 
encounters with the planet during the temporal capture. Using the combined formulae (Equations [30] to 
^ at e = 1, Ae is given by Ci{i)T]^{q)/TQ{q), where TY^{q) and ro(q) are Tk and tq at the pericenter 
distance q, respectively. Therefore, the necessary condition for capture is given by 
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where the interior density of bodies, p^, is as sumed to be 1 g cm the HiU radius of Jupiter is app lied 
to q, and /cap is possibly approximately 100 (iSuetsugu and Ohtsuki Ll2013l : ISuetsugu et al. 11201 ill . As 


shown in Eigurejd] Ci{i) is mainly 0.1 to 10. This density is comparable to or less than the ‘minimum 
mass subnebula’ disk that contains a mass in solids equal to th e mass of current Jovian satellites and 
gas according to the solar composition (ICanup and WardLl2002ll . It should be noted that the temporally 
captured bodies are significantly affected by the central star. However, the temporally captured bodies 
rotate around the host planet, which means that the perturbation by the central star is roughly canceled 
out in a temporally captured orbit. Therefore, the energy loss due to gas drag estimated above may lead 
to bound orbits. 


Inclination decreases during the full capture by gas drag, which is estimated as C' 2 (i) = [{{i)/i)ie/ (e))]e=i 
in EigurejH The initial inclination is damped during capture for 20° < i < 30°, while inclinations re¬ 
main high after capture for other i. 

However, inclinations keep decreasing due to gas drag after capture. A dissipation time of the disk, Tdisk> 
that is shorter than the damping time of inclination is thus necessary for the formation of high-inclination 

























Figure 4 Ci{i) = (Ae)e=i To{q)/TK{q) and C' 2 (i) = ((z)/z(e))e=i derived from the combined equa¬ 
tions (Equations fTdlto [T^. 
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where Mp is the host plan et mass. Since the dissipation processes of circumplanetary disks are not 
clear yet (IFuiii et al. I.I2014I') . it is difficult to discuss the dissipation timescale. However, the dissipation 


timescale needed to form high-inclination satellites seems too short. Therefore, the capture of high- 
inclination satellites might have occurred in the timescale estimated in Equation (1401) before the disk 
dissipation and the resulting satellites tend to have retrograde orbits (see Figure |4]l. 


Summary 

I have investigated the time derivatives of orbital semimajor axis a, eccentricity e, and inclination / of a 
body orbiting in a gaseous disk. 


• I have derived d, e, and z for e > 0.9 and i < H/2a (Equations [T4l to [T^ and for i > 2H/a 
(Equations |30] to |32]l. In addition, I have modified the formulae derived by AHN; Equations (fTTI) 
to (fT3]) are valid for e < 0.2 and i < H/2a, where H is the disk scale height. 

• I have combined the formulae in the limited cases and have constructed approximate formulae for 
a, e, and i (Equations |30]to [38]), which are applicable unless e>lorz>7r — 77/a. 

• Using these formulae, I have discussed the orbital evolution of satellites captured by a circum¬ 
planetary disk. High-inclination satellites are formed if the bodies are captured in approximately 
lO'^ years before the disk dissipation. 
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